Spectral properties of classical two-dimensional clusters. 
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Abstract 

on 



We present a study of the spectral properties like the energy spectrum, 
' the eigenmodes and density of states of a classical finite system of two- 

o ' 

dimensional (2D) charged particles which are confined by a quadratic poten- 



tial. Using the method of Newton optimization we obtain the ground state 



and the metastable states. For a given configuration the eigenvectors and 

a" 

eigenfrequencies for the normal modes are obtained using the Householder 



diagonalization technique for the dynamical matrix whose elements are the 
second derivative of the potential energy. For small clusters the lowest ex- 



^ ■ citation corresponds to an intershell rotation. The energy barrier for such 



rotations is calculated. For large clusters the lowest excitation consists of a 
vortex/anti- vortex pair. The Lindeman melting criterion is used to calculate 
the order-disorder transition temperature for intershell rotation and intershell 
diffusion. The value of the transition temperature at which intershell rotation 
becomes possible depends very much on the configuration of the cluster, i.e. 
the distribution of the particles between the different shells. Magic numbers 
are associated to clusters which are most stable against intershell rotation. 
The specific heat of the cluster is also calculated using the Monte-Carlo tech- 
nique which we compare with an analytical calculation where effects due to 



anharmonicity are incorporated. 
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I. INTRODUCTION 



During the last few years considerable attention has been paid to the study of the prop- 
erties of mesoscopic systems consisting of a finite number of neutral or charged particles. 
The particles are confined by an artificial external confining field. Behavior of either ions in 
a radio-frequency trap (Paul trap) or a Penning trap PJS] and heavy-ion ring storage || can 
serve as an illustration of three-dimensional (3D) Coulomb clusters. Very large Coulomb 
clusters have been created recently in strongly coupled rf dusty plasmas Examples of 
two-dimensional (2D) Coulomb clusters are electrons on the surface of liquid He |5| and 
electrons in quantum dots |§. The vortex clusters in an isotropic superfluid [|7| and in su- 
perconducting grains || have many common features with those of 2D charged particles 
. Refs. JT^JTT|] have been devoted to the investigation of the ground state of 3D clusters 
of charged particles. Below we give a short overview of previous theoretical work on 2D 
clusters of charged particles. 

Clusters of particles in 2D with Coulomb repulsion were investigated by Lozovik and 
co-workers [12| in the case of parabolic confinement. They found that for low temperature 
and in the case of a small number of particles the cluster has a shell structure. A two step 
order-disorder transition was found. With increasing temperature, first intershell rotation 
starts, and intershell diffusion may be possible at high temperature. When the size of the 
cluster is sufficiently large, the simple shell structure gradually disappears in the center and 
features of a Wigner lattice appear. Then cluster melting occurs around the 2D Wigner 
lattice melting temperature. 

Bolton and Rossler [13|] considered the case of parabolic confinement for a small number 
of particles: 1 — 40. They investigated the ground state as well as some metastable states. 
For clusters consisting of 6 particles they determined the barrier height for transition from 
the configuration (1,5) (these are the number of electrons in each shell) to the configuration 
(6). 

Systematic and detailed investigation of the structure of 2D clusters was carried out 



3 



by Bedanov and Peeters ||14|| . They considered both parabolic and hard wall confinement. 



A table of Mendeleev was constructed for clusters with: 2 — 52, 82, 151, 230 number of 
particles. Using the Lindeman melting criterion these authors determined the temperature 
for the order- disorder transition for radial and angular displacement. 

In all of the above works on 2D systems with a finite number of charged particles the 
Monte-Carlo simulation technique was used. We found that in some cases this method is 
rather slow in finding the ground state of the cluster. The reason is that the Monte-Carlo 
technique spends too much time in the vicinity of metastable states such that for a finite 
simulation time not necessarily the correct ground state is found. This becomes more of 
a problem for clusters with larger number of particles which have many more metastable 
states. In Ref. |14| this drawback was partially avoided by heating up the system and cooling 



it down repeatedly. In the present work we will present an alternative approach. To find the 
ground state we choose the Newton method with initial configurations determined randomly. 
In this way we are able to obtain not only the ground state but also the metastable states. 
The latter are relevant in the calculation of thermodynamic properties and the barrier height 
for intershell rotation. 

In previous work flT^-ITl] the ground state properties and melting temperatures were 



obtained. Here we will investigate the spectral properties of the system. This paper is 
organized as follows. In Sec. II we describe the model and introduce the dimensionless 
units. In Sec. Ill our numerical technique to obtain the ground and metastable states is 
outlined and compared to the Monte-Carlo technique. Sec. IV is devoted to the stable 
configurations and the spectrum of normal modes is determined. The barrier height for 
intershell rotation is obtained in Sec. V. Intershell rotation is the lowest excitation for 
small clusters. We correlate the strong dependence of the height of the barrier for intershell 
rotation to the number of particles placed in various shells. In Sect. VI we discuss large 
clusters for which we calculate the density of states and discuss their lowest excitation 
which consists of a vortex/anti- vortex pair. In Sec. VII the zero temperature results for 
the excitation spectrum are used in order to calculate the melting temperatures using the 



Lindeman melting criterion. These results are compared with earlier results |14j] which were 
based on the Monte Carlo simulation technique. As an example of the use of metastable 
states in the calculation of thermodynamic property we calculate the heat capacity in Sec. 
VIII. We compare the Monte-Carlo results with an analytical approach in which we include 
anharmonicity effects in an approximate way. Our conclusions are presented in Sec. IX. 

II. MODEL SYSTEM 



The model system was defined in Ref. 14 . But for completeness we recall the main 



features. Our system is described by the Hamiltonian 

t i>j I '« ' 3 \ i 

where q is the particle charge, e is the dielectric constant of the medium the particles 
are moving in and the confinement potential V(r) = ^mu^r 2 is taken parabolic. Particle 
motion is described by classical mechanics in the plane r — (x,y). To exhibit the scaling 
of the system we introduce the characteristic scales in the problem: ro = (2g 2 /mea;o) 1/ ' 3 
for the length, E = (mco>Qg 4 /2e 2 ) 1//3 for the energy, and T = (mco' 2 g 4 /2e 2 ) 1//3 A;^ 1 for the 
temperature. These scales will be used as our new units and all our results will be given in 
these units. In so doing the Hamiltonian can be written as 



i>j I '» 'i I i 

with V(f) = x 2 + y 2 . The numerical values for the parameters uq, r , E , T for some typical 



experimental systems were given in 14 



In the present paper we will consider only classical systems. Although a classical ap- 
proach for the description of the behavior of electrons in quantum dots is not applicable, 
nevertheless it is possible that certain features of the classical system may survive in a quan- 
tum system. For example in the quantum study of the transition from a crystal to a liquid in 



the absence of a magnetic field ||15|| , we know that the parameter for formation of a Wigner 
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crystal is r s = /o/ a o = 37 ± 5, where Iq is the mean distance between the particles. If the 
number of particles is small, the interparticle distance in the case of parabolic confinement is 
close to ro- Thus for typical parameters for a quantum dot in GaAs with m = 0.067, e = 13, 
hooo = 1 meV we obtain r s = 7.8. Reducing the confinement ujq or applying a magnetic 



field |]TB[ will give us a possibility to investigate the existence of a Wigner crystal or another 



ordered state for a finite number of particles. In Ref. |14| it was found that a classical 2D 



cluster with a finite number of charged particles can be more or less stable than a 2D crystal 
for the same parameter T = q 2 / 'eloksT '. We expect that a similar quantity will be relevant 
in the quantum case and therefore it is expected that also a Wigner crystal like state can 
exist in quantum dots. 

III. NUMERICAL APPROACH 



The Monte-Carlo simulation technique |T7| is relatively simple and provides relatively 
rapid convergence and a reliable estimate of the total energy of the system in cases that 
a relative small number of Metropolis steps is sufficient. However, the accuracy of this 
method in calculating the explicit states may be poor in certain cases. We can understand 
this as follows: for the present system of axial symmetric confinement some configurations 
have very small frequencies for intershell rotation ou m i n = 10 -3 -j- 10~ 4 which may lock the 
simulation in an unstable state. Using the Monte-Carlo method with an acceptable number 
of steps 10 4 -i- 10 5 , in order to limit the computer time, we may obtain the energy E up to 
an error 5, but the error in the coordinates will be proportional to 5 1 ^ 2 /uj rn i n which in such 
a case can be large. 

To circumvent this problem we used a different numerical approach which is mainly based 
on our experience from which we learned that with different modifications to the gradient 
method and the method of molecular dynamics using artificial viscosity we were able to 
obtain more reliable results than with the Monte-Carlo technique. To be more explicit, 
to find the state with the minimal energy we used the modified Newton technique. Since 
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this method is practically not applied in the present field we will give a short outline. Let 
us suppose that the coordinates of the particles in a cluster are given by {r™^ a — x,y, 
i = 1, . . . N} after n-steps in the simulation. Then the potential energy in the vicinity of 
this configuration can be written in the following quadratic form 

H = (r» 4 ) -EE H a,i(r a ,i ~ <J + \ £ £ H a/3jij (r a>i - r^)(r w - rfa), (3) 

i ot ij a,/3 

where H a ^ = —dH/dr a ^ is the force and H a p^j is the dynamical matrix 

d 2 H 

Ha ^ = dr-dr^- (4) 
The next step in our simulation is based on the condition of minimal total energy 

£ XX?7<W>, ij + H a(3 ^{rpj - r n p j ) = H a>il (5) 

3 

where <5 a/ g j y is the unit matrix and the coefficient r\ is added to assure the stability of the 
algorithm. It is easy to show that the iteration procedure converges if rj > — A m .j„, where 
X m in is the minimal eigenvalue of the dynamical matrix. The system of linear equations (|5|) 
is solved using Gaussian elimination. The calculation of the matrix and solving the system of 
linear equations takes about iV 2 numerical operations. This is equivalent to a Monte-Carlo 
step where also about iV 2 operations are needed to find the energy, but the coefficient in 
front of iV 2 is less for the latter. The reason is that to obtain the spectrum of the matrix 
is more laborious. The usual approach guarantees only convergence in the vicinity of the 
minimum. Therefore we introduced an empirical dumping coefficient rj. In the first few 
iterations the value for r\ is set to be large: r\ = 10 -j- 100. If in the next step the total energy 
of the system decreases the dumping coefficient is reduced, while in the opposite case the 
value 7] is increased. From our experience we know that such an algorithm for choosing the 
dumping parameter guarantees convergency of the iteration process. Furthermore, near the 
last steps, the dumping parameter becomes less than the minimal value of the eigenvalue 
of the dynamical matrix and the rate of convergency becomes square ( 5 n+ i ~ <5 2 ). The 
accuracy of the calculated energy 5 is now only limited by rounding errors. For systems with 



axial symmetry there exists an eigenvalue with value zero which corresponds to turning the 
system as a whole around the axis of symmetry. In such a case the second eigenvalue A2 has 
to be taken as the minimal eigenvalue. We found that in order to obtain the configuration 
with minimal energy with an accuracy of | H a ^ \= 10~ 9 -j- 10~ 10 takes about 10 -v- 100 steps, 
the exact number of steps depends on the number of particles. 

After finding the state with the minimal energy we obtain the eigenvalues and eigenvec- 
tors of the dynamical matrix (Q). The eigenfrequencies of it are the eigenvalues squared. 
The condition that the minimal eigenvalue is positive guarantees that the obtained con- 
figuration is stable. Of course also the present method does not guarantee that all stable 
and metastable configurations and the configuration with the lowest energy are found. To 
overcome this difficulty partially we consider a large number (typically 200) of initial con- 
figurations which are generated randomly. From these initial configurations a few stable 
configurations remain, the number of which increases fast when N > 30 -j- 40. Among these 
stable configurations the state with the lowest energy is taken to be the ground state of the 
system. The fact, that usually the state with the minimal energy is achieved already after 
a small number of steps, gives us confidence that this is likely the actual ground state of 
the cluster. Usually, the radius of convergence of the ground state is sufficiently large. We 
confirmed that the present approach for N < 80 leads to the ground state configurations of 
Ref. []14| which were obtained using the Monte-Carlo method with about 10 5 -=-10 6 simulation 
steps. 

The efficiency of the present method is illustrated in Fig.l where we plot the precision 
of the energy, which is defined as the difference from the exact energy value, as function of 
the number of simulation steps for a cluster of 13 particles. It is apparent that the present 
technique converges much faster, about an increase with a factor of 200 is found. Furthemore, 
we discovered that even if within the Monte-Carlo approach the error in the energy is only of 
order 10 -11 , the obtained cluster configuration was unstable. This was found by calculating 
the minimal eigenvalue of the matrix which consists of the second derivative of the potential 
energy with respect to the position coordinates which for the obtained configuration was 
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negative. The present Newton optimization approach did not exhibit such a deficiency. 



In contrast to the Monte-Carlo approach of Bolton and Rossler |13j who found more than 
one stable configuration for the case of iV = 13 particles, the present approach in which 
200 initial configurations were considered, demonstrates that there exist only one stable 
configuration which is (4,9). But for this configuration the minimal excitation frequency 



Umin ~ 6 • 10 is very small which may be the reason for the error in Ref. |13 



IV. EIGENVALUES AND EIGENVECTORS 

A detailed description of the features of the lattice structure, the interparticle distance 
scale in the various shells, and the Mendeleev table for the configurations with N = 2 — 52, 



82, 151, 230 particles was given in Ref. |14[ . Here, we will discuss the excitation spectrum 
corresponding to the ground state configuration of the system. This spectrum is shown in 
Fig. 2 as function of the number of particles for N ranging from 2 to 50. The eigenfrequency 
in this figure is in units of uj J \/2. Notice that there are three eigenfrequencies which are 
independent of N: i) for any axial symmetric system the system as a whole can rotate which 
gives an eigenfrequency uj = 0. This is illustrated in Fig. 3 (figure indicated by k = 1; 
k counts the eigenvalues in increasing order) where the arrows indicate the direction of 
movement of the different particles (i.e. the eigenvectors of the excitation) for a cluster 
with N = 9; ii) there is a twofold degenerate vibration of the center of mass with frequency 
uj = = 1.4142 (see Fig. 3, k — 7 ); and iii) the third eigenfrequency corresponds to a 
vibration of the mean square radius iV -1 J2i( x j + uf) with frequency uj = = 2.4495 (see 
Fig. 3, k = 15). The value of this breathing mode can easily be obtained analytically. 

For clusters of sufficient large size (i.e. N > 8) a typical feature of its spectrum is the 
occurrence of a very low eigenfrequency. Because of the scale in Fig. 2 this frequency is not 
discernable from the uj = frequency and therefore we have listed it in Table I as uj m in- For 
a number of clusters the eigenvectors corresponding to these minimal eigenvalues are shown 
in Figs. 3(/c = 2 ), 4 and 5 (A; = 2 ). For N = 19 and N = 20 the central particle does not 
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move for this specific excitation and consequently its displacement vector has length zero 
and is therefore not visible in Fig. 4. For the clusters with iV = 9, 19 and 20 particles, the 
motion with the minimal frequency w m j„ corresponds to intershell rotation. The necessary 
condition for the existence of intershell rotation is the presence of at least two particles on the 
inner shell in order to conserve total angular momentum. With changing configuration, the 
minimal eigenfrequency can vary by several orders of magnitude (see Table I). For instance, 
for iV = 19 with the ground state configuration (1,6,12), the minimal eigenfrequency is 
u min 0.67, and for iV = 20 and configuration (1, 7, 12), w mill ps 1.0 x 10~ 4 . In both cases 
the minimal eigenfrequencies correspond to intershell rotation (Fig. 4). This large change 
in the size of the minimal eigenfrequency is connected with the shell configuration, and not 
with the total number of particles. For example, if for N = 19 we take the metastable 
configuration (1, 7, 11) whose energy is an amount 1.66 x 10~ 2 larger than the ground state 
energy, we obtain uj m i n w 1.1 x 10~ 4 which coincides practically with io m i n for the cluster 
with 20 particles. 

From the data given in Table I we infer the following law: a high frequency value for 
intershell rotation is obtained for configurations such that the number of particles on the 
outer shell is an integral number times the number of particles on the inner shell. For 
example: iV = 12 (3,9), iV = 15 (5,10), N = 16 (1,5,10), and iV = 19 (1,6,12). For 
clusters with more than two shells (i.e. N > 21) a large io m i n for intershell rotation is found 
for ground state configurations in which the number of particles in the different shells are 
multiples of an integer number. The latter is usually the number of particles in the inner 
shell. For example: iV = 22 (2,8,12), iV = 30 (5,10,15), iV = 45 (3,9,15,18) and to a 
lesser extend also N = 34 (1,6,12,15). These cluster numbers can be considered as the 
magic numbers, because they represent the clusters which are most stable against intershell 
rotations. In previous work by others on 3D clusters magic numbers were determined on the 
basis of energy calculations of the cluster configuration. We found [[14] that for 2D clusters 
no clear steps are found in the cluster energy versus the number of particles in the cluster 
and therefore the stability argument is more appropriate in the present case. On the other 
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hand, a configuration with small u m i n for intershell rotation is realized when the number of 
particles in the different shells have no common denominator. For example: N = 13 (4, 9) 
and N = 20 (1,7,12). 

The above rules can be understood from the Hamiltonian by analyzing the intershell 
interactions using cylindrical coordinates. Let us consider the most simple configuration 
with two shells. From the outset we notice that the occurrence of an particle in the center 
of the system, for example, for a cluster with 20 particles, does not disturb the intershell 
rotation. Therefore, we do not have to consider such a case separately. Let us discuss 
the rotation between two outer shells. The interparticle distance and the distance between 
the shells changes only by a few percent when we alter the number of particles and/or the 
configuration. Therefore, in an initial approximation we can describe each shell by an ideal 
polygon and thus the interaction Hamiltonian between two shells can be reduced to the form 

I Ni N 2 

ff=-E + R l + 2R i R 2 cos (i0i - j9 2 - 6))- 1 ' 2 , (6) 

1 i=i j=\ 

where Ri, R 2 , Q\ = 2ir/Ni, 9 2 = 2ir/N 2 are the radii and angles between particles of the first 
and second shell which have N% and N 2 particles respectively, and 9 is the intershell angular 
distance. The sum (|6|) over the two indexes can be reduced to the sum over one index only 

h = E(^i + R i + 2r i r i c ° s w* - e )y 1/2 > ( ? ) 

where 

0, = (8) 

and I is an integer which is the minimal divider of the number of particles (N\, N 2 , ...) 
in the different shells. From expression ([I]) it follows that the Hamiltonian for intershell 
interaction is periodic in 9. Moreover the period 9*, as a rule is less than the angular 
interparticle distance within a shell. To evaluate the strength of the intershell interaction 
we deduct from the Hamiltonian the following value 
iViiV 2 



Air 



f ^ dx{R\ + R 2 2 + 2RxR 2 cos (x - 9))~ 1/2 , (9) 

J 
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which is independent of 9. This result was obtained from Eq. ([?D by replacing the summa- 
tion over k by an integration. We proved that the error we make in doing so is proportional 
to 91- Numerical summation of ([?[) gives even a more weaker dependence of the interaction 
energy on N for two ideal polygons. When we compare the computed results for the barrier 
height for intershell rotation with those found from Eq.(9) we found that Eq.(9) gives a 
good qualitative description but quantitatively the results are not satisfactory. Therefore 
we may conclude that for small eigenvalues, the exact value of the barrier height is strongly 
influenced by the non-ideality of the polygons. Indeed in order to obtain Eq.(6) we assumed 
that the particles were placed at the edges of an ideal polygon. Because intershell rotation is 
a collective phenomena, one can easily understand that the actual barrier height is less than 
that given by Eq. (|7|) due to the deformation of the polygons during the motion. Indeed, 
during the rotational motion not only the intershell distance changes but also the interpar- 
ticle distance within a shell is altered. This is illustrated in Figs. 3 and 4. From these 
figures we notice that the eigenvectors for the particles in the inner shell have practically the 
same length and are orthogonal to the radius-vector of the particle. For the outer shell the 
situation is different and the eigenvectors have also components in the radial direction and 
futhermore, the length of the eigenvectors are different for the different particles. Therefore 
the vibrations in the radial and axial directions of the outer shell follow the intershell rota- 
tional motion of the particles. Only for clusters in which the number of particle on the inner 
shell is a multiplicative integer factor of those of the outer shell, i.e. when a large intershell 
rotation frequency is found, are the polygons almost ideal which can be understood from 
symmetry reasons and from our numericl results. The characteristics and modelling of the 
intershell rotation will be given in next section. 

When we increase the number of particles, we found that for iV = 39 the motion with 
the minimal eigenfrequency no longer corresponds to intershell rotation, but rather consists 
of rotation of different individual regions of the cluster. For iV > 60 (see Fig. 4) the rotation 
of an inner shell is followed by the rotation of individual polygons at the periphery of the 
cluster. For iV > 115 we found that the minimal frequency io m in no longer corresponds to 
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intershell rotation but corresponds to the excitation of a vortex/anti- vortex pair (see Fig. 5 
for N=151). Higher excitations (see Fig. 5 with k — 4 and k — 6) may consist of multiples of 
such pairs. In case of a cluster of N — 151 particles the 7 th lowest excitation corresponds now 
to an intershell rotation. A more detailed discussion on the nature of low energy excitations 
of large clusters is postponed to Sect. VI. 



V. BARRIER FOR INTERSHELL ROTATION 

In the present section we give a more detailed discussion of the lowest non-trivial exci- 
tation in case of small clusters, which is the intershell rotation. The barrier for intershell 
rotation was obtained using the following procedure. Let us assume that after n-steps in our 
simulation the coordinates of the particles are given by {ff; % — 1, ...N}. After diagonalizing 
the dynamical matrix H a ^^j we obtain the eigenvectors Ai{k) and the eigenvalues = 0Jk- 
The particle coordinates for a next time step is then given by 

2N 

^ +l = ^ + Y.^Mk) ■ (io) 

k=2 

Denote r = r^. as being the lowest frequency for intershell rotation which is taken to be 
constant and which sets the size of the time step. The values of all other coefficients are 
found from the condition of minimal potential energy. This is done as follows: substitute 
the above expression in Eq.(3) which gives us the total energy for the next step 

N 2N oo- -, 2N 

H n+1 = ff„+EEv ^(*) + 9 £ . (") 

i=l k=2 ° r ' 1 A k=2 

from which we readily find the coefficients 

N Q1T 

^ = -VE (12) 

i=l 1 

Trajectories of the particles in four different clusters are depicted in Fig. 6. As was 
mentioned above, intershell rotation takes place in conjuction with radial oscillations. The 
latter are more noticeable for clusters with high symmetry, which have a relative large 
frequency and consequently large barrier heights for intershell rotation. In clusters with 
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only two shells, particles in different shells rotate in opposite direction in order to conserve 
total angular momentum. Such a motion is defined completely by the angle of rotation (p 
of one shell relative to the second. When there are three shells or more it is convenient 
to introduce the angle of rotation *p of the shell with the maximum angular velocity as an 
independent parameter. Fig. 7 illustrates the dependence of the potential energy on this 
parameter ip for two different clusters with N = 9 and N = 40 particles. In general this 
function is well approximated by the simple relation 

,2inp, 



2 



l-cos(— , (13) 

where U*, and (p* are the barrier height and the period, respectively. The values for Z7* and 
ip+ are given in Table I. The above procedure is not able to give the value of the barrier 
height for the clusters N = 39,42,47,51 and 70. The reason is that for those clusters the 
minimal eigenvalue does not correspond to intershell rotation. Notice that the cluster with 
iV = 40 has two minima in the potential energy. One of this minima corresponds to a 
metastable state. Of course U* and ip+ in Table I are determined by the global minimum 
and maximum. In a few cases, the maxima in potential energy are sharper than that given 
by expression ([T3|). The reason is that the energy for clusters with three and more shells are 
not only a function of the angular position of the shell. 

For clusters with two shells, the parameter <p characterizes the motion of the inner shell. 
The angle of rotation of the outer shell relatively to the inner one can be obtained from the 
condition of zero total momentum 

*=(1 + |§V, (14) 

where Ni, Ri are the number of particles and the radius of shell i, respectively. For clusters 
with two shells, the value <p+ presented in Table I is correctly approximated by the simple 
analytical formulas (|8]) and fll4]). The Hamiltonian for intershell rotation taking into account 
kinetic energy can be written in the form 



H = l -N x R\<p\l + §§) + k(l - cos (^)). ! 1.1 ) 

I iV 2 ix 2 2 <p+ 
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For clusters with more than two shells, we can only propose phenomenologic generalization 
to expression fllS]) . Let us label {A^i = 1,...,N} the set of eigenvectors corresponding 
to intershell rotation. Then to first approximation the Hamiltonian for intershell rotation 
becomes 

H= 1 -Rl^+ 1 -U^l-cos( 2 -^)), (16) 
2 2 ip± 

with 

1 M 

fi^TrE^I.l/r, 2 , (17) 

i=l 

where the summation is carried out over the particles of the shell with the maximum an- 
gular velocity. The value of the parameter R± is also given in Table I. Once we have the 
Hamiltonian it is not difficult to find the connection between the barrier angular value (p* 
and the characteristic frequency for intershell rotation 

U* = 2(co min ^^) 2 = 2(u mm S) 2 . (18) 

271 

The parameter 5 has a clear physical meaning: it is the length which a particle travels 
within a shell when it moves over the angle (p*. The approximate expression fllSD is shown 
in Fig. 8 by the solid curve together with the results of our simulation which are given by 
the symbols. Notice that Eq.(|T8"D describes our numerical results very well over an energy 
barrier height variation of more than 8 orders of magnitude. 



VI. DENSITY OF STATES AND VORTEX EXCITATIONS 

From Fig. 2 we notice that the maximum frequency in the excitation spectrum, on the 
average, slowly increases with increasing number of particles. We can easily explain this 
with the aid of the theory of an infinite system. As it follows from our calculations, and has 



been mentioned in previous work |TJj], the minimal interparticle distance decreases slowly 
with the growth of the cluster size due to the compression of the inner shell by particles 
placed at the periphery of the cluster. As a consequence, the maximum value of the wave 
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vector k ~ tt /lo (Iq is the mean distance between the particles) and also the wave frequency 
will increase weakly with the cluster size. 

For large clusters, it is more interesting to consider the density of states (DOS) of exci- 
tations (phonons) which can be obtained by a summation of the energy levels, displayed in 
Fig. 2, over a frequency interval which we took Sou = u max /40, where 0J max is the maximal 
eigenfrequency. The results for iV = 80 and N = 300 particles is shown in Fig. 9. A char- 
acteristic feature in the DOS for all clusters is the occurrence of two broad maxima. From 
earlier investigations [|l8j of classical infinite 2D systems we know that there are two types 
of waves in a 2D Wigner crystal: the lateral sound waves with dispersion relation u ~ k and 
the longitudinal plasma wave with lo ~ y/k, in the long- wavelength limit. Using an analyti- 
cal approximation for the frequency of sound uo\ ~ 0.00363c<jp/i; 2 /Q and the plasma frequency 
lo\ pa u>pklo(l — 0.181483/c/o)? taken from ||18|| , it is possible to show that the positions of the 



broad maxima in Fig. 9 are in qualitative agreement with the ones for an infinite crystal. 
In our dimensionless units u p = 27r/p/o- To obtain the value of u p we used the average 
particle density p = N/rcRl, where R Q is the radius of the most outer shell. The maximum 
frequency of plasma like waves u>2, ma x ~ 1-17oj p for the cluster with A^ = 80 equals 4.67 and 
for A" = 300 is about 5.77. Let us assume that the maximum frequency for sound waves is 
achieved at klo = tt. Then for A^ = 80 we obtain ui iTnax ~ 2.38 and for A^ = 300 we find 
^i^nax ~ 2.94 which are slightly larger than the position of the first maximum appearing in 
Fig. 9. 

From continuum elastic theory, a 2D electron crystal can be considered as incompressible 
at low frequencies ||19|| . In order to check if this is still the case for the present finite 



system we consider the z-component of the rotor ^ r {k) = e z ■ rot^(k) and the divergence 
*&d(k) = divty(k) of the field of eigenvectors of mode k 

1 N 

^d(k) = -^h(k), (19a) 
JV i=i 

1 N 

Mk) = ^ ■ (19b) 

i=l 
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The values ipd,i{k), and ip r ,i(k) for the i particle are given by 



M 

*PdA k ) = E [(fi - r) ■ Mk) + (4 - r) ■ A m (k)] / \ n - r m \ , (20a) 



771=1 



M 



AA k ) = Yl - r) x M k ) + if m - r) x A m (k)\ / | r* - r m \ , (20b) 



771=1 



where r m are the coordinates of the neighbor particles and Ai(k) is the eigenvector of particle 
i for mode k. The rotor and divergence of the eigenvector field are shown in Fig. 10 as 
function of the excitation frequency for clusters of size N = 80 and iV = 300. Notice that 
for small values of the frequency the rotor of the field of eigenvectors is larger than the 
divergence. As a consequence, in a finite system but with iV sufficiently large, the system 
is incompressible and one can expect that the low frequency excitation consists of vortex 
motion in which the particle density is not disturbed. From our computer calculations we 
found that for iV = 151 the minimum eigenfrequency corresponds indeed to the formation 
of a vortex/anti- vortex structure (Fig. 5). Since the total angular moment has to be 
equal to zero, those vortexes always come into pairs. With higher eigenvalues, the number 
of vortexes rises, although this function is not necessarily monotonic (see Fig. 5). Thus 
when iV is sufficiently large the cluster of charged particles can be described as a viscious 
non-compressible fluid in case of small wave vectors. Vortex motion is only expected for 
sufficiently large iV because the velocity of dissipation of the vortex energy is inversely 
proportional to R 2 , where R is the characteristic radius, which increases with increasing N. 



VII. MELTING TEMPERATURES 



In Ref. [inj it was shown that the T = ordered state of the cluster is destroyed with 



increasing temperature (T). The melting temperature for this order- disorder transition 
was obtained by investigating the radial displacement, the relative intrashell and intershell 
angular displacements as function of temperature. Here we will start from the excitation 
spectrum of the zero temperature ordered state in order to calculate the melting temperature 
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using the Lindeman melting criterion |20|. In the harmonic approximation the mean square 



displacement is given by the following expression 



rp 2N 

<^>=T7E^ 2 > (21) 



from which we find the melting temperature T using the relation < u 2 >= 7^, where 
7 = 0.10 for a 2D Wigner crystal E1]J23], and Iq is the mean interparticle distance. As 



discussed in previous section, there exists a number of configurations with uj min very small 
which will give a very large contribution to the sum (pTj) . In order to see what the effect 
is of these very low frequency excitations on the melting temperature we also considered 
the sum without the first term. Then we will find the temperature for intershell diffusion. 
Because for large clusters, the value of the interparticle distance around the center and near 
the periphery can be considerably different, therefore we will use the mean value of relative 
displacement in order to define the melting temperature 

N 2N 



T = 7A 



i=l k=2 



-1 



(22) 



where Ai(k) is the displacement vector for the i-particle in mode k, and U is the mean 
interparticle distance for the ith particle. 

The numerical results are shown in Fig. 11. As we expect there is a significant difference 
in the transition temperature whether intershell rotation is taken into account or not. These 
results agree qualitatively with the results of Ref. |14| where it was found that: i) for clusters 



with a small number of particles the angular order is destroyed at much lower temperatures 
than the radial order which agrees with the large difference in melting temperatures shown in 
Fig. 11; ii) for larger clusters both temperatures are practically equal as is also apparent in 
Fig. 11. Orientational order and radial order disappear practically at the same temperature 
for A > 40; and iii) the melting temperature at which intershell diffusion sets in, is a 
decreasing function of the number of particles in the cluster up to about A ~ 20 -j- 40 
beyond which it starts to increase, which agrees qualitatively with Fig. 11. The magnitude 
of the transition temperature found in the present approach is slightly higher than found 



in the Monte Carlo study of Ref. ||14|| . This is a consequence of the present harmonic 
approximation which has a limited validity near the melting temperature. The melting 
temperature for intershell rotation (top part of Fig. 11) is strongly influenced by the value 
of uJmin which is proportional to the rigidity of the cluster agains intershell rotations. In fact 
it is a measure of the stability of the cluster against intershell rotations. As was mentioned 
before the value of w m j„, and also U+, is determined by the configuration of the cluster. 
Clusters with a magic number of particles have a large melting temperature for intershell 
rotation. These fine details were not present in Ref. flLH| . 

It is known that for an infinite crystal the sum (|2~2"D diverges logarithmically in the low- 
wavelength which is due to the presence of lateral sound waves. Therefore one uses the 
average square displacement of interparticle distance in Lindeman's melting criterion. In 
our case, such criterion gives the relation 



T = jN 



N 2N i M 1 _1 

E'r a E^ a iE(4(*)-A»(*)) a 

Li=l k=2 m=l 



(23) 



where the sum over m runs over the M neighbor particles. The numerical results obtained 



using Eq. (|23|) is shown in Fig. 12. These results are very close to those found in Ref. [[bf] with 
the exception that here near N ~ 150 a maximum is found while the Monte Carlo results 
slowly increases towards the N — > oo value. We want to emphasize that if the number of 
particles is not too large, the transition temperature obtained with the second criterion fl23|) 
is lower than the one from Eq. (|22]) . This indicates that the particles mainly move towards 
each other, and only for N > 200, the effect of small wave vectors begin to appear. In the 
latter case the neighbor particles move with the same velocity and the difference in the value 
of critical temperatures obtained using the spectrum of the eigen vibrations (Fig. 11) and 
the Monte-Carlo technique |TJ] is very small. 



VIII. SPECIFIC HEAT 



Before we already mentioned, that in order to obtain the ground state we generated 
many initial configurations and in so doing not only stable states but also metastable states 
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were obtained. Thus, if we also calculate the spectrum of normal modes u>k for each local 
minimum we can easily obtain the partition function within the harmonic approximation 
and consequently all the thermodynamic quantities like the free energy, the specific heat,... 
Such an approach was followed in Refs. |23|] and |24]] for 3D clusters where also the influence 
of anharmonicity and saddle points on the partition function was studied. In Refs. PSipj 



only, the main characteristics of the spectrum of normal modes was used. Here we know the 
complete spectrum of our finite 2D system and are therefore able to calculate the partition 
function more correctly. 

In the quasiclassical approximation the partition function for a cluster with N particles 
is given by |f25| 



Z{T) = (2nh)- 2N J dqdpexp (— H (q, p) / k b T) , (24) 

where q = (r 1; r N ), p = (pi,...,p N ) are 2N dimensional vectors. The partition function 
can be written as 

M 

Z(T) = £ exp {-U m /T)Z m {T) 7 (25) 

m=l 

where Z m is the partition function of the mth metastable state whose energy differs with 
the ground state by an amount U m . The dimensionless units for temperature and specific 
heat are used here and below. In the vicinity of this mth metastable state, the Hamiltonian 
is quadratic in the normal coordinates. Because the energy barrier for intershell rotation 
is small, the effect of anharmonicity will already appear at low temperatures. Therefore 



we will integrate only over a small region of particle motion | qi |< yj2U m (k)/T(Vk, m which 
results in 

27V y 

Z m = 9mZ rot J] 7 erf{^U rn {k)/T), (26) 

where Z rot oc VT is the part of the partition function resulting from the rotational degrees of 
freedom, g m = 2^/6+ is the degeneracy of the mth state, which is determined by the number 
of particles occupying a shell, U m (k) is the barrier height for normal mode k, and erf(x) 
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is the error function. These parameters are given in Table II for a number of metastable 
states. 

For convenience let us consider only one normal mode. At low temperature T <C U m (k) 
expression (p6| ) results in the usual value for the specific heat for a harmonic oscillator 
C = 1. For high temperature T <C U m (k) the specific heat equals 1/2 as for free motion. 
For the intermediate temperature region T ~ U m (k) expression ( ^q) gives an interpolation 
between these two limiting cases. Unfortunately, we know only the value of the barrier 
height for intershell rotation. For the remaining normal modes, we will use the analogy with 
the Lindeman criterion to write the phenomenological relation 

U m (k) = luN 1 ^ , (27) 

where Iq is the mean interparticle distance, and 7„ = 0.2 4- 0.3. The above expression is 
then used in the numerical evaluation of the partition function ( p6|) and (|25|). Below we will 
mainly deal with the specific heat 

d 2 d\ogZ 

C -&t T ~&T' (28) 

which is shown by the solid curve in Fig. 13. 

General features of the behavior of the specific heat as function of T and iV can be 
predicted without detailed information regarding the metastable states and the spectrum 
of the normal modes. At low temperatures such that T <C (U m (k), U m ) the specific heat 
is only determined by the ground state, and the effect of anharmonicity is not essential. 
Consequently C = 2N — 1/2 as is apparent in Fig. 13.. Usually the barrier for intershell 
rotation is the smallest energy, which is also smaller than the difference in energy between 
the ground state and the metastable states. In the temperature range U^fc = 2) C T < 
(Ui(k 7^ 2), U m ) the specific heat will be constant and having the value C = 2N — 1. 
Such a small reduction in C is visible in Fig. 13 near T ~ 10~ 3 . With further increase of 
the temperature, the behavior of the specific heat is determined by the competition of two 
processes. On the one hand, transitions to metastable states which lead to an increase of 
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the specific heat, and on the other hand, the effect of anharmonicity which will reduce C. 
This interplay will lead to peaks in the specific heat as is apparent in Fig. 13. Note that 
the position of the peak does not equal the melting temperature. 

In the order- disorder transition region the applicability of the above approach is ques- 
tionable. Therefore, we also calculated the specific heat using the standard Monte-Carlo 
technique. As the initial state we took the ground state of our system. Then we fix the tem- 
perature and execute 10 5 steps of the Metropolis algorithm to allow the system to achieve 
equilibrium. Next about (4 4- 10) • 10 6 steps of the Metropolis algorithm are made in order 
to reduce the statistical error. The specific heat is then found using the following formula 

C = N+{<E 2 >-<E> 2 )/T 2 , (29) 

where E is the potential energy for the system with N particles. In Fig. 13 we compare 
the results from the Monte-Carlo simulation (full dots) with the above results (full curve) 
which are based on the excitation spectrum of the T = stable and metastable states. 
Note that for the small cluster with N = 9 very good agreement is obtained. For the other 
two clusters good quantitative agreement is found at low temperature while at intermediate 
and high temperatures only the qualitative behavior is correctly described. Thus for large 
clusters the approximate model is not able to give a satisfactory description of the effect of 
anharmonicity. Nevertheless there is qualitative agreement in the position of the maxima. 
We have tried to vary the parameter 7„ and to change the integration interval for allowed 
particle motion in Eq.(26) but we were not able to obtain any better agreement. 

IX. CONCLUSION 

We have presented the results of a numerical simulation of the ground state and the spec- 
trum of normal modes of classical 2D clusters with quadratic confinement. The barriers for 
intershell rotation and the specific heat are also obtained. The Lindeman melting criterion 
in conjunction with the T = excitation spectrum of the ground state configuration was 
used to obtain the order- disorder transition temperatures for angular and radial melting. 
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For systems with axial symmetry, and an intermediate number of particles the normal 
mode with the lowest frequency corresponds to intershell rotation if there are at least two 
shells. A low excitation energy for intershell rotation is found for clusters which have a shell 
configuration such that the number of particles on each shell have no common multiple. If 
the number of particles in the outer shell is an integer multiple of the number of particles 
in the inner shell, the cluster will be most stable against intershell rotation which define 
the clusters with magic numbers. Such clusters also have a large melting temperature for 
intershell rotation. Distortion of the axial symmetry of the external potential, will lead to a 
rise in the eigenfrequency and in the barrier height for intershell rotation. For large clusters, 
i.e. N > 100, the normal mode with the lowest frequency corresponds to a vortex/anti- 
vortex excitation. 
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FIGURES 

FIG. 1. Accuracy of the calculated groundstate energy versus the number of simulation steps 
using the Monte-Carlo technique and the present optimized version of the Newton technique for a 
cluster consisting of N = 13 particles. 

FIG. 2. Excitation spectrum of normal modes as function of the number of particles in the 
cluster. The frequency is in units of uj/y/2. 

FIG. 3. Eigenvectors for a cluster with N = 9 particles for different mode number k which 
correspond to the eigenfrequencies u\ = 0, ll>2 ~ 0.127, ujj « 1.414, W15 ps 2.449. 

FIG. 4. Eigenvectors corresponding to the minimal frequency for clusters with N = 19, 20, 39, 
and 60 particles. 

FIG. 5. Eigenvectors for the cluster with N = 151 particles for four different values of the mode 
number k. 

FIG. 6. Trajectories of the particles for intershell rotational motion in case of clusters of 
N = 9, 19, 20 and 38 particles. 

FIG. 7. Cluster energy versus shell turn angle for the shell with the maximal angular velocity 
for clusters with N = 9, and 40 particles. 

FIG. 8. Energy barrier of intershell rotation versus (uj5) 2 where u> = uj m i n is the frequency for 
intershell rotation and S is the linear distance traveled by a particle over one angular period. 



FIG. 9. Density of phonon states for clusters with N = 80 and 300 particles. 

FIG. 10. The value of the divergence (a,c) and rotor (b,d) of the displacement field for clusters 
with N = 80 (a,b) and N = 300 (c,d) particles. 

FIG. 11. Temperature of cluster melting obtained using the Lindeman criterion with and with- 
out intershell rotation. 
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FIG. 12. Melting temperature for large clusters as obtained from Lindeman's melting criterion 
excluding intershell rotation and incorporating the relative displacement of neighbor particles. 

FIG. 13. Specific heat for clusters with N = 9, 25 and 34 particles as function of temperature. 
Solid lines are the results as obtained from an approximate calculation of the partition function 
and the solid dots are results from the Monte-Carlo simulation. 
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TABLE I 



Table I. Shell configuration (N 1 , N 2 , . . .) for clusters with N-particles with parabolic con- 
finement. The minimal excitation frequency (to m i n in units of u /^/2), the period (</?*) in 
degree units and the barrier height ([/*) for intershell rotation are given together with the 
parameter R± for the ground state of the cluster. 
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TABLE II 



Table II. Shell configuration (A^, N 2 , . . .) for some metastable states for a number of 
different clusters. U m is the energy difference of the metastable configuration with the 
ground state energy and W m = J7 Wk,m=i/ U ^k,m is the relative statistical weight. 
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